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1. Introduction 

Dual Scattering Chamiel (DSC) schemes are characterized by a two-step 
cycle of iteration which alternately updates the computed fields within cells 
and on their interfaces. If the updating instructions are explicit, then a 
near-field interaction principle leads to the typical structure of the DSC 
algorithm and to a scattering process interpretation. A pair of vectors that 
represent the same field within a cell and on its surface essentially consti- 
tutes a scattering channel. Equivalently, scattering channels are sometimes 
defined as pairs of distributions that 'measure' the field within the cell and 
on one of its faces. A well known DSC scheme is the Transmission Line Ma- 
trix (TLM) numerical method along Johns' line |JoB| wherein transmission 
line links visualise the scattering channels. The TLM picture of wave prop- 
agation fails in Computational Fluid Dynamics, thus giving rise to the far 
more general DSC setup. 

In technical depth, DSC schemes and their relation to the TLM method 
have been analysed in [Helj . Another study shows that they are uncondi- 
tionally stable under quite general circumstances, made tangible with the 
notion of a-passivity |He2j . DSC schemes are particularly suitable for han- 
dling boundary conditions, non-orthogonal mesh, or also for replacing a 
staggered grid where otherwise need of such is. 

In section |2 we first recapitulate some characteristic features of DSC 
schemes, which in extcnso are treated in |Hel| . before wc sketch in section|21 
the Oberbeck-Boussinesq approximation to the Navier-Stokes equa- 
tions. The DSC model then outlined in sectionQlrepresents a prototype, first 
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of all. In fact, the Boussinesq equations for viscous quasi-incompressible 
flow, inspite of retaining the often predominant non-linear advcctive term of 
the Navier-Stokes momentum equations, can only claim limited range of 
validity, due to their well known simplifications. The Oberbeck-Boussinesq 
approach is yet prototypical also in providing a basis for many turbulence 
models |ATP| which can be implemented following the general guide lines 
of this paper. 

The algorithm with non-orthogonal hexahedral mesh cell outlined in sec- 
tion^has only recently been implemented by Spinner, where it is directly 
coupled to a Maxwell field solver, allowing thus for computing conductive 
and convectivc heat transfer simultaneously with the heat sources of a lossy 
Maxwell field. Besides the underlying ideas and simplifications that enter 
the Oberbeck-Boussinesq approximation and its DSC formulation, a set 
of numerical results are displayed to illustrate and validate the approach. 

2. Elements of DSC schemes 

Given a mesh cell, a port is a vector valued distribution associated to a cell 
face which assigns a state vector = {p , Z) to a physical field Z (any 
suitable smooth vector valued function in space-time). 

We also require that, within the given cell, a nodal image p~ of p exists, 
such that 

(1) Z) = [po a, Z) = {p,Z o a-'), 

for every Z (of class C°°, e.g.) , where a denotes the spatial transla- 
tion cr : — > R"^ that shifts the geometrical node (centre of cell) onto the 
(centre of) the respective face, cf. Figl. 




Fig. 1. Port on a cell face with nodal image. 

DSC fields split thus into port and node components, and z" , which 
represent the field on the cell surfaces and within the cells. The two compo- 
nents are updated at even and odd integer multiples, respectively, of half a 
timestcp r and arc usually constantly continued as step functions over the 
subsequent time intervals of length r 

Moreover, we assume that the updating instructions arc explicit, i.e., 
with possibly time dependent functions F and G , for t = mr ; m G N 
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where [z]t stands for the entire sequence up to time t 
[z]t := {z{t - ^T))^gN 

and we agree upon fixing zP'"(t) := for t < 0. (The 'back in time 
running' form of the sequence has certain technical advantages, cf. |Helj .') 

A fimdamental DSC principle is near-field interaction, which is related to 
the Courant-Levi stability criterion in spelling that every updated state 
depends only on states (up to present time t ) in the immediate neighbour- 
hood. More precisely: The next nodal state depends only on states (along 
with their history) in the same cell and on its boundary, and a subsequent 
port state depends only on states (with history) on the same face and in 
the adjacent nodes. 

As a consequence of near-field interaction, every DSC process allows for 
an interpretation as a multiple scattering process in the following sense. 

If M is a mesh cell system and d( denotes the boundary of cell C G M, 
then every DSC state permits a unique scattering channel representation in 
the space 

with canonical projections tt^' " : P —> P^' " into the port and node com- 
ponents of cell C ( the cell index can usually be omitted without danger of 
confusion). Also, there is a natural involutary isomorphism nb : P —> P 

nb : {zP, zP^ ) ^ (zP", zP) , 

which is called the node-boundary map and obviously maps PP onto P" 
and vice versa. For every DSC process z = (zP , z"'){t) , the following 
incident and outgoing fields zf^^ and z"„j are then recursively well defined, 
and are processes in and P" , respectively: 

For t < , zf„ (t) : = 2:"„t(i— §) : = , and for every < t = mr ; 
m e N 

zPit) - n6ozl,(t - I) , 

z-{t+^)- nbozlit)- 

At every instant holds, hence, zP (t) ~ n6 o z"^^^ ( t — ^) + zP^{t) and 
z" (t -t- ^ ) — nb o zP^{t) -\- (t ^ ) . Then, near-field interaction 
implies that every state is only a function of states incident (up to present 
time t ) on scattering channels connected to the respective node or face. 
More precisely, by induction holds 

Theorem 1 . A pair of functions H and C exists, such that for every cell 
C, ^ M the process z ? = tt p o z complies with 



(3) 



(0 
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and the port process = o z satisfies 

+ r) ^ e((^li(^ + ^ -M^))«|9C;MeN) • 

(5) ^ 

( ' I ' short-hand for 'adjacent to' ) 

Remarks 

(i) The statements immediately imply that z " and z^ -^^ are themselves 
functions of states incident on connected scattering channels, since 

^C%«t(*+o) = 3?((^m(^ -M^))p6aC;/^eN) - and 

(6) ^ r ~ r 

(ii) 51 and 6 arc named the reflection and connection maps, respectively, 
of the DSC algorithm. 

(iii) Near field interaction implies computational stability, if the reflection 
and connection maps are contractive or a-passive |He2| . in addition. 



3. The dynamic equations 

The DSC algorithm is thus simply characterized as a two-step explicit 
scheme that alternately updates states in ports and nodes of a cellular mesh 
and which, in virtue of a near-field interaction principle, allows for a canon- 
ical interpretation as a scattering process, the latter exchanging incident 
and reflected quantities between cells and their interfaces. 

Ports and nodes are related to physical fields by vector valued distri- 
butions that evaluate the fields on cell faces and within cells of a cellular 
mesh. Such a distribution may be a finite integral - as is the case with the 
TLM method, where finite path integrals over electric and magnetic fields 
are evaluated in a discrete approximation to Maxwell's integral equations 
|He8| . In the simplest case, it is a Dirac measure that pointwise evaluates 
a field ( or a field component ) within a cell and on its surface. The dis- 
tribution can also be a composite of Dirac measures that evaluate a field 
at different points in the cell - which applies, for instance, to the gradient 
functional treated in sect.^J 

Classical thermodynamics with, in particular, energy conservation entail 
the convection- diffusion equation for the temperature T in a fluid of veloc- 
ity u with constant thermal diffusivity a , heat source(s) q , and negligible 
viscous heat dissipation, viz. 

dT 

(7) h u • grad T = aAT + q . 

at 

This is the energy equation for Boussinesq-incompressihle fluids, e.g. jGDN| . 
The Navier-Stokes momentum equations for a fluid of dynamic viscosity 11 , 
under pressure p , and in a gravitational field of acceleration g require 
d 

(8) — ( gi u ) + ( u ■ grad ){qu) + gradp = /i A u + g . 
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The Oberbeck-Boussinesq approximation |Obbj . [5ssj starts from the as- 
sumption that the fluid properties are constant, except fluid density, which 
only in the gravitational term varies linearly with temperature; and that 
viscous dissipation can be neglected. Equations © become so with Qoo = 
const and g (T) = g oo P {T {t — Too ); P ■= g^^ d g / dT 

(9) ^ + {u-grad)u + = A u + /3 ( T ( x ) - Too ) g . 

ot goo goo 

Applying the Gauss-Ostrogradski theorem to the integrals over A ~ div grad 
on cell C with boundary and using u • grad f + f div u div ( / u ) . 
equations (|3 ^ yield, with a time increment r , the following updating in- 
structions for nodal T and u, the latter, along with q , averaged over the 
cell volume V"^ 

T{t+^) : = T + T{Tdivxi + q) + 

(10) T r 

+ — / {agradT - T u) ■ dF 
^< Jdc 

and 

u(i + J) : = u + Tiudivu - /3(T - r^)g - ) + 

^ Qoo 

(11) ^ r I. 

+ —{ — gradu ■ dF - u ( u ■ dF ) } 

JdC Qoo JdC 

In equations H1(J|I and Hll|l the last former updates (at time t — r/2) of the 
nodal quantities enter the right-hand sides in the first line, and the former 
updates (at time t) of the cell face quantities enter the second line. Note 
that divu vanishes, of course, for incompressible flow. 

All nodal values of T and u are thus updated at the reflection step 
of the algorithm, while the cell face quantities that enter the integrals at 
the right-hand sides are updated on the connection step. In next section we 
will see how to proceed with an unstructured (non-orthogonal) hexahedral 
mesh. 



4. Non-orthogonal hexahedral cell 

The physical interpretation of a DSC algorithm associates a smoothly vary- 
ing ( e.g. in time and space C°°- ) scalar or vector field Z to port and node 
states and of a mesh cell system. 

Let any hexahedral cell be given by its eight vertices. Define then edge 
vectors (j/e)i/=o,...,ii, node vectors (^6)/i=o,i,2, and face vectors (t/)t=o,...,5, 
using the labelling scheme of figure |21a 

1 x-^S 

— T > M ^ ,e M = 0,1,2 

(12) (-1)' 

^ ' and J := ( ,,+,„e +,,^,,^^(_^,,,, e) A 

^ ((4+20^ +(6+20 e) i = 0,...,5, 
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with all indices understood cyclic modulo 12 and A denoting the cross 
product in R^. At every cell face l G {0, 5} and for any given t G K_|_ the 
following time shifted finite differences of Z in directions ^5 {fj, = 0, 1,2) 
form a vector valued function 

R , . i 2(-lY(Z'' u_^,2- ,ZP It) ifu=[t/2l 

( [ x ] denotes the integer part of a; G R ) . The time increments are chosen 
conform with the updating conventions of DSC schemes (as will be seen in a 
moment) and are consistent. In fact, in the first order of the time increment 
r and of the linear cell extension, the vector Z in the centre point 
of face 6 approximates the scalar products of the node vectors with the 
gradient VZ . Let, precisely, for a fixed centre point on face l and e G R+ 
the e-scaled cell have edge vectors ^e'^ : = e . Let also Z^ denote 

function l(T^ for the e-scaled cell (with node vectors = C/^^ )• Then at 
the fixed point holds 

(14) < ^6,grad(Z) > = ,,b-\/Z = lim lim 1,V^~Z^, 

as immediately follows from the required C^-smoothncss of the field Z. 

To recover the gradient VZ from l|13() in the same order of approxima- 
tion, observe that for every orthonormal basis (i/u),y=o....,m-i of R™ or C™ , 
and for any basis (^j&)^i=o,...,7Ti-i with coordinate matrix = < yU , i^fi >, 
the scalar products of every vector a with equal 

^ — ^m — 1 

(15) < ^b, a > = 2^ _p < ^b , ,yU > < yU , a > ^ jS'^a^, 

( at the right-hand side, and henceforth, we observe Einstein's summation 
convention - without yet summing up over indices that anywhere appear as 
left-hand subscripts ) , hence 

(16) a. =7^^ with ili:):^m)T' ■ 

I.e., the scalar products of any vector with the basis vectors ^5 transform 
into the coordinates of that vector with respect to an orthonormal basis 
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by multiplication with matrix 7 = , where = < jyii , > , i.e. 

(3 is the matrix of the coordinate (column) vectors with respect to the 
given ON-basis nU , and 7 its adjoint inverse. 

This applied to the node vector basis and p4|l yields the approximate 
gradient of Z at face l 

(17) ,VZ, = 7^; ,V^Z^. 

The scalar product of the gradient with face vector i^f^ = < ,/, yU > , 
V e {0, 1,2} is thus 

= : tS*' 

Continuity of the gradient at cell interfaces yields linear updating equations 
for ZP on the two adjacent faces. In fact, for any two neighbouring cells C, 
X with common face, labelled u in cell C, and k in x , continuity requires 

(19) [S = -Is. 

Substituting H18() for ,''5' and ^5 and observing the time shifts in H13|l 
provides the updating relations for ZP at the cell interfaces. To make these 
explicit, we first introduce the quantities i,zP' " , ( t = 0, 5 ; /i = 0, 1, 2 ) 



(20) Xit) 



2{-iy z" \t = [l/2] 

{2^+iZP - 2^^ZP)\t-r/2 else 



which in virtue of yields ,,z^ = {p, Z) = {p^, Z o i,a^^ ) = 

= i,z^^ I Z o , where : n ^ p denotes the nodal shift pertinent to 

face i . In particular 

(21) .zlf,^{t) = 2(-l)\ZP|,, 
which together with (|20|) is consistent for ^ ^ [t/2] with 

(22) ,z;(i + r/2) = -\{2^+iZ^^ + 2^.zjl){t). 
From (1121 Ha El) follows that 

^S\t+r = .s''(.<|*W2-2(-l)'4'/'l.^^*+.) 

= .s^(,z;:i,+./2-<5['/2i,zP|,+.). 

The continuity of Z , i.e. \zp ^ ^ZP , the implies with H19I20|) 



c p ^ :g^:^;^(t + r/2) ^^:^(t + T/2) 

For completeness we agree upon setting ^zf^it + t) := \z'^{t + t /2) 
for yu 7^ [t/2] (which yet contains a slight inconsistency, in that continuity 
might be infringed; this can be remedied by taking the arithmetic means of 
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the two adjacent values). - In fact, our agreement doesn't do harm, since 
any discontinuity disappears with mesh refinement. 

We have, hence, a complete set of recurrence relations for { given z " 
by the former reflection step) which at the same time determine the field 
components on face t and their gradients 

(25) ,VZ, = 7^;,zP. 

Essentially this constitutes the connection step of the algorithm. 
Nodal gradients are similarly (yet even more simply) derived using 

V^^;(t+^) := (2^+1^^ - 2M^P)(t); M = 0, 1, 2 

in the place of (|13|l and then again (|17|1 . With the node and cell-boundary 
values and gradients of T and u the nodal updating relations for the latter 
are immediately extracted from equations f llOllllI ') in sect. 13. For equation 
this is essentially (up to the convective term) carried out in |Hel| . 
sect. 5, and the procedure remains straightforward in the case at hand. 
Note that a well-timed LES coarsening routine |He3j should be periodically 
carried out before the nodal step of iteration in order to avoid instablities 
from the energy cascade \Po\ - 

The updating relations thus obtained are explicit and consistent with 
near-field interaction (only adjacent quantities enter). So, they can op- 
tionally be transformed into scattering relations for incident and reflected 
quantities Q along the guidelines of section|21- with established advantages 
for stability estimates [He2| . 

5. Pressure 

Conservation of mass simply requires divergence-free flow, div u = , for 
a Boussinesq-incompressible fluid. In integral form, using Gauss' Theorem, 
this means = div u dV — Jg^ u • dF . Since equations ([7| |H1) in sec- 
tion|31do not a priori guarantee this, additional arrangements must be made. 
This is done by the following procedure, which in Magneto-Hydrodynamics 
is known as divergence cleaning. 

In a successive overrelaxation (SOR) routine, carried out between the 
connection and reflection steps of the iteration cycle, firstly the (discrete) 
cell boundary integrals /g^ = J^^u-dF are computed and then the pressure 
p which compensates lec so that 

(26) — / gradp ■ dF = / u ■ dF . 

Taking indeed Z of the preceeding section as the pressure, equations H26|l 
( of course, in discrete form with sums over the cell faces ) yield a unique 
solution p" = for every cell, given the right-hand side integral Idc, ■ 
Note that we are actually solving Poisson's equation Ap = ( Qoc/t) div u 
in integral form. 

With the new face pressure gradient computed in following the lines of 
section 0] the face values of u are updated as u — {t/ goo) gradp . 
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After each SOR eycle, continuity of grad p at the cell faces is restored 
by updating the port values of p according to the instructions of the last 
section. The loop of processes is reiterated until J2q ^9C ^ ^ ^ suitable 
bound e (which happens after a few iterations for appropriate choices). 

6. Convection in coaxial line 

We illustrate the approach in a stalwart application by displaying the re- 
sults of simulations with coaxial line RLlOO-230 under high power operating 
conditions ( which are realistically inferred from a ion cyclotron resonance 
heating ICRH experiment in plasma physics). 

The inner and outer conductors of diameters 100 mm and 230 mm are 
made of copper and aluminium, respectively, and the rigid line is filled 
with air at atmospheric pressure. The heating process has been simulated 
from standby to steady state CW operation, at frequency 100 MHz and 160 
kW transmitted power, for horizontal position of the line and with outer 
conductor cooled at 40 degrees Celsius. 

Figure 13 b displays the computed air flow profile (vertical section) in 
steady state, which is attained some minutes after power-on. Visibly, the 
natural convection pattern is nicely developed. 

Our computations have been carried out with a 3D-mesh of 10 layers in 
axial direction, over 200 millimeters of line, the transverse cross section of 
which is displayed in figure 13 a. At the metallic interfaces no-slip boundary 
conditions arc implemented and free-slip conditions at all other boundaries. 

Simultaneously, and using the same mesh cell system, a Maxwell field 
TLM algorithm was run to provide the heat sources. 




7. Conclusion 

A prototypical implementation of the Oberbeck-Boussinesq approxima- 
tion to viscous flow has been presented in this paper, which demonstrates 
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the fundamental fitness of the DSC approach for fluid dynamic computa- 
tions. DSC schemes thus significantly transcend the range of application of 
the TLM method from which they descend. 

A next natural step in the line of this study is the implementation of turbu- 
lence models which are compatible with the Boussinesq approach, such as 
the k — e model jATPj , first of all. Another objective is the incorporation 
of compressible flow. 

- We hope this paper stimulates some interest into joint further investigation 
in these directions. 
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